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This report documents the recent developments in methodologies for the evaluation of the in- 
tegrity and durability of composite structures, including i) the establishment of a stress-intensity- 
factor based fracture criterion for bimaterial interfacial cracks in anisotropic materials (see Sec. 2); 
ii) the development of a virtual crack closure integral method for the evaluation of the mixed-mode 
stress intensity factors for a bimaterial interfacial crack (see Sec. 3). Analytical and numerical 
results show that the proposed fracture criterion is a better fracture criterion than the total energy 
release rate criterion in the characterization of the bimaterial interfacial cracks. The proposed vir- 
tual crack closure integral method is an efficient and accurate numerical method for the evaluation 
of mixed-mode stress intensity factors. 


1 Introduction 

Laminated composites made of graphite epoxy have a higher stiffness modulus than aluminum, 
as well as other prefered material properties, which make them highly attractive. As a result, 
their application have been steadily increasing in the past decades. Among the various failure 
mechanisms for laminated composites, the propagation of a delamination crack is often the leading 
cause of global failure for a structural component. Therefore, the potential disbonding is one 
of the critical problems in the evaluation of the integrity of composite structures. This research 
is to enhance the capability to predict the strength of structural components made of composite 
materials. 

The prediction of the onset of propagation of a delamination crack in a laminated composite 
was often based on the total energy release rate of the delamination crack. However, various 
recent experiments have pointed out that the critical energy release rate of a delamination crack is 
dependent on the combination of individual modes of failure; mode I (opening), mode II (shearing), 
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and mode m (tearing). For example, the critical energy release rate for a mode I fracture is an 
order of magnitude less than the critical energy release rate for a mode III fracture in a T300- 
5208 unidirectional laminate. Hence, the critical total energy release rate for any two delamination 
cracks would be quite different unless the ratio of mixity of the fracture modes is similar. This 
indicates that total enegy release rate is not a good candidate as a fracture criterion for a mixed 
mode fracture. 

Some researchers have suggested the use of the components of mixed-mode energy release 
rates as the fracture criteria for a delamination crack between two dissimilar lamina. However, 
numerical (finite element) results [Sun and Jih (1987) and O’Brien (1982)] indicates that the in- 
dividual components, in mode I and mode II, of the energy release rates for an interfacial crack 
are strongly dependent on the mesh refinement of the finite element model. Raju et al. (1988) has 
shown mathematically that the mixed-mode energy release rates of an interfacial crack does not 
converge upon mesh refinement, using crack-closure integrals in the context of an analysis based 
on finite element method. In addition, Chow and Atluri (1995) have shown that the three individual 
mixed-mode energy release rates, ( G/, G//, G/// ), are insufficient to characterize the stresses and 
displacements near the interfacial crack tip. Coupled mixed-mode energy release rates, ( G//_///, 
Gj-m, Gi-n ) are required to fully describe the asymptotic stresses and displacements near an in- 
terfacial crack tip. Therefore, it is cumbsome to use components of energy release rates for fracture 
criteria. 

However, the mixed-mode stress intensity factors of an interfacial crack are not dependent 
on the numerical model, unlike the mixed-mode energy release rate components. In the area of 
interfacial fracture mechanics for isotropic bimaterials, the mixed-mode stress intensity factors 
are preferred by many researchers as the fracture parameters to define the toughness function of 
an interfacial crack. The mixed-mode stress intensity factors for an interfacial crack between 
dissimilar isotropic media were first defined by Williams (1959), and Rice and Sih (1965). They 
have found that the stress oscillates and the crack surface overlaps near the interfacial crack tip in a 
linear elastic fracture analysis. However, this solution can be justified [Rice (1988)] in engineering 
sense when the nonlinear material behavior and crack surface contact zone exists only in a small 
scale at the crack tip. 

The acceptance of stress intensity factors for laminated composites has been slow, until recently 
a clear mathematical definition of stress intensity factors for an interfacial crack in an anisotropic 
bimaterial continuum has been established [see Wu (1989); Qu and Li (1991) and Qu and Bas- 
sani (1993)]. This definition of stress intensity factors reduces to the classical definition of stress 
intensity factors when the bimaterial continuum degenerates to a homogeneous one. This report 
documents the development of a fracture criterion based on this definition of stress intensity factors 
and numerical methods for its evaluation. 
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2 A Stress-Intensity-Factor Based Fracture Criterion 
2.1 Crack tip solution for a bimaterial interfacial crack 

Consider a two dimensional deformation in an anistropic material, in which the displacements 
depend on the Cartesian coordinates, x\ and X2 only. The stress cr ij and strain £/./ satisfy the 
following constitute relation. 


3 3 

j = ^ Q ( ^ ) 
*= 1/=1 

Three distinct complex eigenvalues pj (with positive imaginary parts), and two eigenmatrices 
A 3 x 3 and B 3^3 can be defined from the modulus’ Qju using Eq. 2. 
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Qik — Qlfcl Hik — Cj 1 k2 

The eigenmatrices A and B can be normalized so that 


Tik = Caki 


( A7 BT ] [ .1 1 0 3 ] [ B ] = l6 - 


( 2 ) 
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where I is the identity matrix. 

Two real matrices, S and L, can then be defined by the normalized eigenmatrices, A and B, as 

L = — 2iBB r S = i (2AB r - 1) (5) 

In an anisotropic bimaterial medium of material ( 1 ) and material (2), a bimaterial matrix f) can be 
defined as 

P= +l( 2 ) -1 ) -S (2) L (2)_1 ) (6) 

where superscripts (1) and (2) signify material #1 and material #2, respectively. 

Using the complex variables, Zj (j=l,2,3), defined by 


Zj=x\+pjX 2 (7) 

one can express a displacement field (u,) and its assoicated stress field (o ji) in the anisotropic 
medium using three analytic complex functions, f\ (zi),/ 2 ( 22 ) 5 / 3 ( 23 )- 
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Uj = 2Re 
CT], = — 2Re 
G 2 ( = 2Re 


[ fj{tj ) 

L;=l 7 
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U =1 


I^fe) 

L / =1 


( 8 ) 

(9) 

( 10 ) 


Consider a semi-infinite and traction-free crack between two anisotropic media of material (1) 
and material (2) [see Beom and Atluri (1994) for more details]. The three complex functions 
fif\zj), for material #1 and #2 are: 


where 


/ '5 2,< ' ,fe) = ^^ B<2r '( I -^ Y( ^ )8 " 

Y(zf) =I+2p { z f ~ z i“j ^ + j)2 [* _ 2 { P 

e= 5“ ln (r?f) 


( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


and &j (j=l,2) are the normalized base vectors. Note the difference between the matrix P and the 
scalar variable p. 

When the analytic functions, fj ; (zj), are substituted into Eq. 8-10, one obtains the crack tip 
stress field and displacement field. 

The stress along the interface is 

%{x ' )= zm v(xi ' )k 

where T is { (J 12 , ^22>^23 } r ; and k = {Kn,Ki,Km} T is the vector of stress intensity factors. 


( 16 ) 
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2.2 Equivalent stress intensity factors for the characterization of mixed-mode singular stress 

fields 


Consider an interfacial crack between two laminae, [0,a], where the ply angles of the upper and 
lower plies are 0 and a respectively. For a crack in a homogeneous material, in which 0 is equal 
to a, the singular stress field along the interface can be decoupled into three individual modes: <722 
governed by Kj\ <5\2 governed by Ku\ and <723 governed by Km. 

However, such one-to-one relation does not exist for a bimaterial interfacial crack. As shown 
in Eq. 16, the oscillation index, e, produces all the stress components along the interface for each 
mode governed by K ), Ku or Km. Therefore, the singular stress field cannot be decoupled into the 
three unique individual modes. Each of the singular stress fields along the interface is a function 
of the three stress intensity factors coupled by the oscillation index. 


It is postulated that the fracture of an interfacial crack is determined by the singular stress field 
in the region near critical damage radii. The equivalent stress intensity factors fc ro = {£//, i£/, £///} 
at a characteristic radius r a can be defined to characterize the singular stress field in the region near 
critical damage radii. 1 

t r „=Y(ri e ) k (17) 


Then, 


T(*l) = 


1 Y(*f) l_ fi 

V27aTY(/£) ~v/2twT 


(18) 


It is seen from Eq. 18 that the equivalent stress intensity factor £ makes the coupling very weak 
near the characteristic radius r 0 . 

The damage radius for a mode I and mode II crack can be approximated with 



Using the experimental data from Tab. 1, the damage radius of an interfacial crack in a T300- 
5208 graphite-epoxy laminate are ri = 0.05 mm and r// = 0.24mm. Since the magnitude of these 
two damage radius are close to the magnitude of the ply thickness (t p i y = 0. 127mm), the ply thick- 
ness, tpiy, can be used as the characteristic length used in the defintion of the equivalent stress 
intensity factor, fc. 

With the ply thickness as the characteristic length, the coupling between (<712 and £/) and 
(<722 and £//) for the interface [0/90] are significantly reduced near the failure region. Along the 


1 It can be shown that the definition of (i r „ is invariant of different measuring units [Rice (1988)]. 
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Table 1 : Experimental data for the critical energy release rate and their e- 
quivalent stress intensity factor of T300-5208 



G/c 

G//c 

(J/m2) 

Gmc 

Kic 

Knc 

(MPa y/m) 

Kmc 

°22 °12 
(MPa) 

Ref. 4 

103 

279-587 

- 

0.91 

2.7-3.8 

- 

- 

- 

Ref. 5 

88 

154 

- 

0.84 

2.0 

- 

- 

- 

Ref. 6* 

- 

- 

1200 

- 

- 

2.3 

- 

- 

Ref. 22 

- 

- 

- 

- 

- 

- 

50 

75 

Average 

96 

340 

1200 

0.88 

2.9 

2.3 

50 

75 


* Experimental value for AS4/3501-6 graphite epoxy. 


Stress Singularity Near A [0/90J Interfadal Crack 



ct b due -In 
<i a doe to£ « 1 

a a due to * 1 

a a due to it, *! 



Figure 1 : The decoupling of stress near the interfacial crack using ply thickness as the characteristic 
length for stress intensity factor definition 
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interface, 0.1 t p i y < x\ < 10 t p i y , the normalized coupling between the singular stress field and the 
stress intensity factor is less than 10% as shown in Fig. 1. As a result of this decoupling effect, the 
equivalent stress intensity factors of an interfacial crack has a similar physical interpretation with 
the stress intensity factor of a homogeneous crack near the damage region, in the sense that the 
strength of the singular stress field along the interface near the damage radius for mode I, mode II, 
and mode in are governed by Rj, Ru, and Rjh respectively. 


<*22 C*i) ~ 

Ol2(*l) ~ 

<723 (*l) « 


\/2nx\' 

Rn . 

y/2KXi ’ 

Riu_ 

\/2nx\ 


for 


X\ — tply\ ^ ^ 0 


( 20 ) 


2.3 A proposed stress intensity factor based fracture criterion 

Rice (1988) first suggested that the interface fracture toughness be defined by the total energy re- 
lease rate and the phase angle of the stress intensity factor, \|f = tan -1 {Kjif Kj). Since then, many of 
the fracture toughness experiments for an interfacial crack are based on these two fracture param- 
eters. Wang and Suo (1990) have used the Brazil nut specimen to obtain the interface toughness 
function for a Plexiglas/epoxy interface. Shanbhag et al. (1993) performed some fracture tests for 
perspex/epoxy, aluminum/epoxy and steel/epoxy interfaces using the compact tension specimen. 
Yuuki et al. (1994) have used the Brazil nut specimen to obtain the interface toughness function 
for an aluminum/epoxy interface. However, instead of using the energy release rate and the phase 
angle as the fracture criterion, Yuuki suggested a fracture criterion based on the quadratic mixture 
of the normalized stress intensity factor: 



= const 


( 21 ) 


where K\c is the critical value for a pure Ri and Rue is the critical value for a purel£}/. By 
normalizing the individual components of the stress intensity factors, the relative importance of 
each component is physically more intuitive. Unlike the fracture criterion based on the phase 
angle, this criterion can be easily extended to include the component Rm without the loss of the 
physical insight. 


In the present work, the quadratic mixture of the stress intensity factor is modified to include 
the component Rjjj. Furthermore, an additional factor of 0.85 for Kjc is added to reflect on the 
experimental observation [Yuuki et al. (1994)] that K IC drops significantly when a small Ku/Kj 
exists as shown in Fig. 2. With this additional factor, the quadratic mixture theory would be valid 
for a wide range of Kn/Ki and this theory would deviate from the experimental results only when 
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Figure 2: Fracture results for an interface crack 


Ku/Ki is very small. The fracture criterion used to predict the onset of a delamination crack growth 
is as follows: 



£o > 1 — > interfacial crack failure 


( 22 ) 


2.4 Interfacial fracture criterion for laminated composites 

Consider an example of a delamination crack between the [0/90] laminate as shown in Fig. 3. The 
0 lamina denotes that the fiber direction is along the out-of-plane direction, while the 90 lamina 
denotes that the fiber direction lies parallel to the crack. Since all of the experimental fracture data 
for the critical stress intensity factor, Kc, is based on a delamination crack in a unidirectional lam- 
inate, a question arises on whether the fracture data for [0/0] or [90/90] laminate could be used to 
predict the onset of delamination crack in the [0/90] laminate. Fig. 4a shows the stress distribution 
along the interface between the [0/90] laminate near the crack tip; and Fig. 4b shows the stress 
singularity near a crack in a homogeneous laminate ([0/0] or [90/90]) with the stress intensity fac- 
tor oiki/Kic = 1. These figures show that the stress distributions for both the cracks (bimaterial 
and homogeneous) are very similar except for a small region where there is considerable material 
nonlinearity and damage (defined as Zone 1 in the figures). Since the fracture of the delamination 
crack is assumed to be based on the stress concentration near the damage radius and the stress field 
near this damage radius for the interfacial crack is very similar to the homogeneous crack, the crit- 
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Interfacial crack between [0/90] laminate 



>x. 


Figure 3: The critical stress intensity factor for an edge delamination crack between (a) [0/90] 
laminate, (b) [90/90] laminate, and (c) [0/0]laminate 
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ical stress intensity factor for [0/90] laminate, £[ 0 / 90 ]c< should be the related to the value of either 
^[o/o]c or ^[90/90]C- In addition, Lucas (1992) has performed some experiments on the fracture of 
a delamination crack in unidirectional laminates of different ply angle and found that the fracture 
toughness of a [90/90] laminate is less than a [0/0] laminate, ^[ 9 o/ 90 ]c < ^[ 0 / 0 ]c- Since a delami- 
nation crack between two dissimilar lamina can fracture at either the upper or lower lamina, it is 
reasonable to assume that the delamination crack would fracture at the material of lower fracture 
toughness. Hence, using this assumption, the critical stress intensity factor for [0/90] laminate, 
£[ 0 / 9 o]c> would be postulated as^[ 9 o/9o]c- 

The nature of the stress field near a bimaterial interfacial crack is often defined by mixed-mode 
stress intensity factors. To deal with this, Chow and Atluri (1996) have suggested a mixed-mode 
fracture criterion based on the quadratic mixture of the mixed-mode stress intensity factors: 

(ofc) +(fe) + (fc) =( 0 /«A) 2 (23) 

£ 0 >£c^ interfacial crack failure 

where Qf ac is a non-dimensional constant dependent on the angle between the crack front and 
the fiber’s direction. From the experiments performed by Lucas (1992), Qf ac for [90/90] laminate is 
1.00, Qf a c for [45/45] is 1.15, and 0/ ac for [0/0] laminate is 1.61. This experimental result indicates 
that the [90/90] laminate has the lowest fracture toughness followed by the [45/45] laminate and 
then the [0/0] laminate. Using the assumptions discussed above, in which the bimaterial interfacial 
crack would fail at the material with a lower fracture toughness, the constant Q/ ac for [45/90] and 
[45/0] would have the same values as [90/90] and [45/45] respectively. Therefore, Qj ac for [45/90] 
and [45/0] would have the values of 1.00 and 1.15 respectively. These two values would be used 
in the post-buckling analysis discussed in the next section. Since the normalized stress intensity 
factor, J^o, is composed of mixed-mode stress intensity factors normalized with their respective 
critical stress intensity factors, ko is a non-dimensional value. In the examples of composite lam- 
inates under uniaxial tension, Chow and Atluri (1996) have demonstrated that the critical stress 
intensity factor, kc, has the value of one when the delamination crack propagates in a 2-D fashion 
in which crack front remain straight. However, in the post-buckling failure of a stiffened laminat- 
ed composite panel, the delamination crack propagates in an elliptical shape rather than in a 2-D 
fashion. As a result, the critical stress intensity factor, kc, to be used in the post-buckling analysis 
would have to be determined from the experimental data. 


2.5 The effect of characteristic length 

Chow and Atluri (1995) have studied the effect of characteristic length on the normalized stress 
intensity factor, ko . In the few examples they have considered, the value of ko changes by only a 
few percent when the characteristic length changes by an order of magnitude. The relationship of 
mixed-mode stress intensity factors for different characteristic lengths is given by Qu and Bassani 
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a singular stress field in a [0/90] laminate with Ki/Kjj^ 901 = 1 
(ft, is defined with r 0 = tpiy) 



b singular stress field in a homogeneous laminate 
( [0/0] or [90/90] ) with K|/k(J )/901 = 1 


Figure 4: Stress distribution near the (a) interfacial and (b) homogeneous crack tip with Zone 
1 being the damage zone ( or nonlinear region) and Zone 2 being the region where the stress 
concentration influences the onset of crack propagation. Since the stress concentration at zone 2 
are the same for both laminates, the fracture criterion for both laminates should be similar 


(1993) as 


fe ri =Y(r,/ro)fe ro (24) 

/S C /S /V /V T 

where Y is a bimaterial matrix for the interfacial crack and k = < Ku,Kj,Kjh j 

Consider a delamination crack between two angle-ply lamina, q and f, and the stress intensity 
factors, k ro , are varied in a unit sphere: 
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<f> 


Figure 5: A contour plot of maximum change in the value of for various ply lay-ups [0, <()], when 
the stress intensity factor, k rQ , are varied in the unit sphere with the characteristic length reduced 
by an order of a magnitude 



(^olr 0= 0.1tp|y ^°lro“lp*y^ 0 l r O-tp*y 

Figure 6: The probability of any combination sets of stress intensity factors, {ku,fci,kuj} T of a 
delamination crack in any [0/<J>] laminate which would result in the value of ko to be changed by a 
certain factor when the definition of the characteristic length is reduced by an order of a magnitude 
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fci = sin (a) cos (P) 

fcii = sin (a) sin (P) (25) 

km = cos (a) 

where - 180° < a, P < 180°. By varying the stress intensity factors, k ro , in the unit sphere defined 
in Eq.25, a contour plot of the maximum change in the normalized stress intensity factor, k 0 , for 
various ply layups -90° < 0,<j> < 90°, when the choice of characteristic length is reduced by an 
order of a magnitude is shown in Fig. 5. Except for the region near [0/90] interface, the effect 
of the characteristic length is quite small. Even for the [0/90] interface, the maximum change in 
k 0 is less than 13%. To put things in perspective, a few percent change in ko is really not very 
significant from an engineering point of view given the fact that most of the experimental data 
for delamination fracture would vary by 10%-20%. Moreover, the number of combinations of 
(k//,k/,km} T and [0/<(>] which gives significant change in the value of ko is quite small; for 90% 
of the possible combinations of (k//,k/,km} T and [0/<J)], the value of k 0 would change by no 
more 5% as shown in Fig. 6. Only 0.5% of the possible combinations would cause the value of k 0 
to differ by more them 10%. 


2.6 Validation 

The analysis presented in this section is based on the quasi-three-dimensional finite element method 
under generalized plane strain condition. This finite element analysis is based on eight-noded 
quadrilateral, isoparametric elements with 3 degrees of freedom per node. Standard quarter-point 
elements are used to model the interfacial crack tip. The numerical calculations performed in this 
section are for the T300-5208 graphite epoxy using the following elastic properties: 

Ei = 137 GPa, Ej = Ez = 10.8G.Pa, Gzt — 3.36 GPa, Gzl = Gjl = 5.65 GPa, Vjz = 0.49, 
Vzz, = vrz, = 0.238, t piy = 0. 127mm 

Here, L is the longitudinal direction, T is the transverse direction, and Z is the direction through 
the thickness. The experimental data for the critical energy release rate and the stress intensity fac- 
tors of a mode I, II and HI crack are listed in Tab. 1 . These experimental results were obtained from 
double cantilever beam specimens with an initial crack front normal to the fibers’ direction. In the 
present section, the critical stress intensity factors are based on the average of these experimental 
data. The fracture criterion used to predict the onset of delamination crack growth is based on the 
quadratic mixture of the normalized stress intensity factor defined in Eq. 22. 


Mode I and II dominant crack 

The stress intensity factor for a delamination crack in a laminate under tensile loading can be 
mode I and II dominant when the stacking sequence consists of 90° plies being laid-up with plies of 
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Table 2: Laminates with mode I and II dominant delamination crack 


Laminate 

Crack between 

£33 £/ £33 £// 

Kic Kiic 

£ 33^7// 

Kmc 

total 
G IC 

[±30/ ±30/90/90] 5 

-30/90 

213 86 

-11 

5.8E4 

[±45/0/90] 5 

0/90 

146 29 

-1.4 

1.8E4 

[±452/02/902] s 

O 2 / 9 O 2 

206 41 

-2 

3.5E4 

[± 453 / 03 / 903 ] s 

O 3 / 9 O 3 

253 50 

-2.4 

5.4E4 

[0/ ±45/90] s 

-45/90 

80 65 

-6.9 

1.8E4 

Laminate 

Crack between 

Predicted 

Real 

Error 

Laminate 

Crack between 

Failure §33 

Failure £33 

£33 -£33 
§33 

[±30/ ±30/90/90] 5 

-30/90 0.0038 

0.0035 

8 % 


[±45/0/90] 5 

0/90 

0.0057 

0.0056 

3% 

[± 452 / 02 / 902 ] s 

O 2 / 9 O 2 

0.0041 

0.0047 

- 12 % 

[±45 3 /0 3 /90 3 ] 5 

O 3 / 9 O 3 

0.0033 

0.0037 

- 10 % 

[0/ ± 45/90] s 

-45/90 

0.0087 

0.0080 

3% 


different angles. The mismatch of the lateral contraction due to the Poisson’s effect between the 90° 
plies and plies of different angles can sometimes be sufficient to cause the onset of delamination 
crack growth. Some of the examples for mode I and II dominant crack are listed in Tab. 2. 

Tab. 3 lists the calculated stress intensity factor and energy release rate for two free-edge delam- 
ination cracks in a [±30/ ± 30/90/90] s laminate. The criterion based on the total energy release 
rate predicts that the crack would propagate along the second [+30/ - 30] interface. However, the 
criterion based on stress intensity factor predicts the delamination crack would propagate along 
[—30/90] interface which agrees with the experimental observation from O’Brien’s research. 

In Tab. 2, the criterion based on the total energy release rate predicts both [0/ ±45/90] s and 
[±45 /0/90] s laminates would fail under the same tensile load since both laminates have a delam- 
ination crack of the same energy release rate. However, the criterion based on the stress intensity 
factor predicts the [0/ ± 45/90] 5 laminate would fail at a strain 50% greater than the [±45/0/90] 5 
laminate which agrees with the experimental results. Furthermore, Tab. 2 shows that the predicted 
failure strains, based on the stress intensity factors, of various laminates are within 12 % from the 
experimental results [O’Brien (1982) and Sendeckyi et al. (1975)]. Hence, it validates the use of 
stress intensity factor as the fracture parameter to predict the onset of delamination growth. 

Note that all the delamination cracks listed in Tab. 2 lie on top of a 90° ply in which the normal 
of the crack front is in the fibers’ direction. Moreover, the normal of the delamination crack front 
on the double cantilever beam specimen, in which the critical stress intensity factor is obtained, 
is also in the fibers’ direction. Therefore, the critical stress intensity factor listed in Tab. 1 can 
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Table 3: Delamination cracks on a [±30/ ± 30/90/90] s laminate 


Interface crack between 

E33 £/ 

K,c 

£33 kji 

Kuc 

£ 33 ^/// 

Kmc 

£0 

^33 Gfotal 
Gic 

second +30/ — 30 

143 

55 

108 

207 

7.0E4 

-30/90 

213 

86 

-11 

265 

5.8E4 


be directly used in the examples listed in Tab. 2. However, the critical stress intensity factors do 
change if the normal of the crack front deviates from the fibers’ direction. The experiment by Lucas 
(1992) has shown that as the angle between the normal of the crack front and the fibers’ direction 
increased, the critical stress intensity factor would increase as well. Hence, some adjustment on 
the fracture criterion would be required when the interfacial crack does not lie either on top or 
bottom of a 90° ply. 


Mode III dominant crack 

When a laminate of ±0 layers, where 0 is between 0° and 90°, is under tensile loading, the coupling 
between the normal strain, 633, and the transverse shear stress, CT23, are directly opposite for +0 
and —0 layers. As a result, the stress intensity factor of a delamination crack between ±0 layers 
is often mode m dominant when the laminate is under tensile loading. In this study, the fiber 
orientations of \0°, 30°, and 45° are considered. The stacking sequence of the laminate is either 
[(±0) 2 ] s (alternating) or [©2/ — ©2] s (clustered). 

Tab. 4 lists the calculated stress intensity factors and the comparisons between the predicted 
failure strain and the experimental results by Herakovich (1982). There are several explanations 
that can be given for the considerable disparity between these results. For [(± 10) 2 ] s and [(±30) 2 ] s 
laminate, Herakovich observed no delamination propagation when these laminates failed. In fact, 
using the Tsai-Hill failure theory, we can see that these laminates failed because the stresses on 
each lamina exceeded the ultimate stress. The reasons for the prediction errors for other laminates 
can be attributed to the wrong assumption used in describing the crack front. Observations from 
the experiments for these stacking sequences show that the delamination propagates in a triangular 
shape rather than in a straight crack front shown in Fig. 8. Hence, to fully model the delamination 
growth for these laminates, a full three-dimensional model must be used instead. Furthermore, 
materia] nonlinearity must also be taken into account in analyzing the failure of ±45° laminate. 
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Table 4: Laminates with mode m dominant delamination crack 


Laminate 

Interface crack between 

£33 %i £33 £// 
Kjc Kiic 

£33^7// 

Kinc 

ZllGtotal 

Qjc 

[(±10) 2 ] 5 

second ±10/ — 10 

-2.6 2.2 

100 

3.9E4 

[(±30) 2 ] s 

second ±30/ — 30 

11 12 

87 

2.3E4 

[(±45) 2] s 

second ±45 / — 45 

6.9 8 

31 

3.8E4 

[± (IO2) ] s 

±10 2 /-10 2 

-8 0.4 

173 

1.2E4 

[± (30 2 ) ] s 

±30 2 /-30 2 

-7.5 0.0 

160 

7.3E4 

[± (45 2 ) ] s 

±45 2 /-45 2 

-2.6 0.0 

61 

1.3E4 

Laminate 

Interface crack between 

Predicted 

Real 

Error 



Failure §33 Failure £33 

£33 -£33 

§33 

[(±1°) 2] 5 

second ±10/— 10 

0.0100 

0.0070 

43% 

[(±30) 2] s 

second ±30/ — 30 

0.0113 

0.0117 

-3% 

[(±45) 2 ] 5 

second ±45 / — 45 

0.0301 

0.0158 

90% 

[± (10 2 ) ] s 

±10 2 /-10 2 

0.0058 

0.0049 

18% 

[± (30 2 ) ] s 

±30 2 /-30 2 

0.0063 

0.0067 

-7% 

[± (45 2 ) ] s 

±45 2 /-45 2 

0.0163 

0.0097 

68% 


3 A Virtual Crack Closure Integral Method for the Evaluation of Mixed-mode SIF for 

Interfacial Cracks 

This section documents the development of a virtual crack closure integral method for the evalua- 
tion of mixed-mode stress intensity factors for an interfacial crack between two dissimilar isotropic 
and anisotropic mediums. The mixed-modes stress intensity factors can be related to a set of en- 
ergy release rates by the use of a coupled mixed-mode energy release rate. Therefore, the stress 
intensity factors can be easily solved for once the mixed-mode energy release rates are evaluated 
from a finite element solution. 

The modified crack closure integral method has become a useful technique in calculating the 
fracture parameters of a crack lying in a homogeneous isotropic domain. This method is attrac- 
tive for many researchers and engineers because of its simplicity and effectiveness in calculating 
the individual components of the mixed-mode stress intensity factors and the energy release rates. 
The Irwin’s virtual crack closure method was first applied by Rybicki and Kanninen (1977) with 
4-noded isoparametric finite elements to compute the mixed-mode energy release rates of a crack 
lying in a homogeneous isotropic domain. To improve on the accuracy of the computed frac- 
ture parameters, Raju (1987) has further extended this method for higher order and singular finite 
elements. 

To calculate the stress intensity factors for an interfacial crack in an isotropic bimaterial con- 
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tinuum, Lin and Mar (1976) have developed the hybrid stress element method, and Kathiresan and 
Atluri (1978) have developed the hybrid displacement element method. Using the mutual integral, 
M, Cho et al. (1994) have calculated the stress intensity factors for an interfacial crack between t- 
wo dissimilar isotropic media from the stress results of an ordinary finite element analysis (without 
the need to develop a hybrid element). Chow et al. (1995) have further extended both the hybrid 
stress element method and the mutual integral method to compute the stress intensity factors for 
an interfacial crack between dissimilar anisotropic media. However, these methods required the 
knowledge of the asymptotic stress and displacement field around the interfacial crack tip. Since 
the asymptotic stress and displacement field involve complex numbers, the calculation of the stress 
intensity factors using these methods require a considerable programming effort. 

The computation of the stress intensity factors for an interfacial crack can be significantly 
simplified, using the virtual crack closure integral approach. The stress intensity factors for an 
interfacial crack can be computed from the crack opening displacements and the nodal forces 
at and ahead of the crack tip. The simple formula to calculate the stress intensity factors for 
an interfacial crack are given in three separate categories: an isotropic bimaterial continuum, an 
orthotropic bimaterial continuum, and an anisotropic bimaterial continuum. 


3.1 Virtual crack closure integral method 

The work required to extend a crack by an infinitesimal distance A is equal to the work required to 
close the crack to its original length [Irwin (1958)]. Thus, the energy release rate for mode I and 
mode II deformations can be expressed as 

G, = ~ / A a 22 (r) 62 (A - r) dr (26) 

2A Jo 

GlI= TX f a i2 ( r )$i (A — r)dr 
2A Jo 

where <7 is the stress distribution ahead of the crack tip, and 5 is the crack opening displacement 
behind the crack tip. Here, the crack-axis is along the x\ axis and normal to the *2 axis as indicated 
in Fig. 7. 

The procedure to obtain these energy release rates from finite element solutions were given 
by Rybicki and Kanninen (1977) and Raju (1987). For quarter-point singular elements (shown in 
Fig. 7), the energy release rate for mode I can be obtained from the nodal forces and displacements 
near the crack tip (along the crack-axis) with the virtual crack closure technique as 



( 27 ) 
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*2 


m 


Figure 7: Quarter-point singular elements surrounding the crack tip 



where 

3k 

tu = 6 — — ; tn = 6tc — 20; (29) 

1 

*21 = 2 ’ t22 = i ’ 

The symbol represents the nodal force computed from the singular elements in material 1 in 
the direction parallel to X 2 axis at node [i] as shown in Fig. 7. 

For a crack in a homogeneous domain, there exists no coupling between {022 (r) and 5i (A — r) } 
and between (ci 2 (r) and 52(A-r)}. However, these couplings do exist for an interfacial crack 
between two dissimilar media. Because of such couplings, the calculated energy release rates G/ 
and G// are dependent upon the value of the distance A. As a result, Raju et al. (1988) as well as Sun 
and Jih (1987) have found that the calculated energy release rates using the virtual crack closure 
integral method are dependent upon the mesh refinement of the crack tip elements surrounding 
an interfacial crack. To resolve this problem, Chow and Atluri (1995) has introduced a coupled 
energy release rate G/_// to simplify the relationships between the stress intensity factors and the 
mixed-mode energy release rates (obtained through virtual crack integral method). The coupled 
energy release rate G/_// is defined as: 

G/_//=— [ [<7i2(0 $2(A — r )+<*22(f)Si (A — r)] dr (30) 

2A Jo 

Gi -11 can be evealuted for the nodal forces and displacments at the nodes of the quater-point 
singular elements at the crack tip. 
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Once these energy release rates are calculated, the stress intensity factors for an interfacial 
crack can be easily computed from a linear equation. The detailed relations for between stress in- 
tensity factors and energy release rates for (i) an isotropic bimaterial continuum, (ii) an orthotropic 
bimaterial continuum, and (iii) an anisotropic bimaterial continuum, are discussed in the remaining 
part of this section. 


3.2 Interfacial crack in an isotropic bimaterial continuum 

For an interfacial crack between dissimilar isotropic media, there exists an oscillation index which 
causes the oscillation of stress and displacement near the crack tip. This oscillation index is defined 
as 


£ = 



/g) +/iWk( 1 > 

ifi) +^(0k( 2 ) 


(32) 


where p is shear modulus and k = 3 - 4v for plane strain and k = (3 - v)/ ( 1 + v) for plane stress. 
The stress field along the interface at a distance r ahead of a crack tip is given by 


a 2 2 + ;a ]2 = -f=r , ' e (33) 

V2nr 

where K = K/ + iKn is the complex stress intensity factor. The crack opening displacements at a 
distance r behind the crack tip are given by: 


62 + i'Si = 


cW+cP) 

2\/2k ( 1 + lie) cosh (tee) 


K^/?r i£ 


where the compliance parameters, c, are 


c 0) = 


k( j ) + 1 
/|0> ’ 


cP) = 


k( 2 ) -1- 1 

—fFT 


(34) 


(35) 
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Substituting Eq.33 and Eq.34 into Eq. 26-30, the energy release rates can be defined as func- 
tions of the complex stress intensity factor, K. 


c(’)+c( 2 ) 

Gi = — — r lim/?e 


87t cosh (tie) 


A 2iZ K 2 


:h + 


\K[ 


1 -j- 2/e 1 -I - 2/e 


-h 


(36) 


c 0) +c (2) 

Gu = — c lim/?e 

87tcosh (7te) 

c (0 -p c ( 2 ) 


A 2iz K 2 r 
nzh + 


\K\‘ 


1 + 2/e 1 -I - 2/e 


-h 


Gi-n 


47tcosh (tie) 


lim Im 


A 2i£ K 2 T 
—I l 


1+2/e 


(37) 

(38) 


where 


I\ = / cos 2 co(sin(ocosco) 2 ' e </o) (39) 

Jo 

h = J^ 2 cos 2 0) (cos 0)/ sin G)) 2 **/© = (40) 

The integral I\ in Eq.39 is a Beta function and can be evaluated numerically using numerical 
integration scheme. To simplify the relations between the stress intensity factors and the energy 
release rates, the integral I\ is approximated with 


71 


cosh (Tie) 


(i-k-W 


(41) 


This approximation is valid for the realistic range of e. Dundurs (1969) has shown that for an 
interfacial crack in isotropic bimaterial continuum, the absolute value for the oscillation index, e, 
has to less than or equal to ( 1 /2 tt) In 3. The error of the approximation in Eq. 41 for this range 
normalized with J^ 2 cos 2 G)//£D is less than 0.5%. 

Substituting Eq. 39-41 into Eq. 36-38 and rearranging the equations, the simple relationship 
between the stress intensity factors and the energy release rates are obtained: 


K 2 1 ai\Gi — a\iGi-u 1 

Kjj ) \ a2\Gn + a\2Gi-u ) 

where 


a 11021 + fl 12 fl 22<321 ~ a h 

«22«21— «12 <211^21 +«12 


]{ 


(42) 


32 cosh 2 (tic) (l +4e 2 ) 
cO) +c( 2 ) 


( 43 ) 
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an = (1+4 e 2 ) + ^1 - cos[2eln(A)]+ ^y£- sin[2eln(A)] (44) 

«12 = 2 j ^yE-^E 3 ^ cos [2e In (A)] - ^1 sin[2eln(A)] J (45) 

fl22= (1+4e 2 )- ^l-^E 2 ^cos[2£ln(A)]- ^y£- jE 3 ^ sin [2 e I n (A)] (46) 

« 2 i = 4 j^l -^£ 2 ^cos[2£ln(A)]+ ^yE-^£ 3 ^ sin[2£ln(A)]| (47) 

Hence, the square of the stress intensity factors, Kj and Kj h can be easily solved from the 
linear equation in Eq.42. The signs for the stress intensity factors, Kj and Kjj, can be ascertained 
from the crack opening displacements. For a crack in a homogeneous isotropic continuum, E = 0, 
resulting in a considerable simplification of Eqs. 43-47, i.e.: 

Kj = - G /, Kjj = -G n (48) 

c c 

3.3 Interfacial crack in an orthotropic bimaterial continuum 

For an orthotropic bimaterial aligned with the interface crack coordinate system (e.g., a 0/90 degree 
laminate where the fibers are parallel to the crack front in one material and perpendicular in the 
other), the stress field along the interface at a distance r ahead of a crack tip is given by Qu and 
Bassani (1993) as: 


„./ 2 „ , m i« p! / 2 *,+;p! / 2 a:„ « 

P 2 022 + .P, 0,2- ~ j = r 

and the crack opening displacements at a distance r behind the crack tip are 



s 2 +/p; /2 5 l 


(s^ — s^ 2 )) \/2nrr Jt 
TtP ( 1 + 2 /e) cosh (tie) 


(p ^Ki + i^Kn) 


The bimaterial constants 


(49) 


(50) 


E = 



(51) 


P = \/Pi P2 


(52) 
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P.= 


$(') —s( 2 ) 

1 y(l)p(')/'T|(0 -f" j(2)p(2) /t|(2) 


(53) 


P 2 = 


c(l)_,(2) 


.?(') /rj ( 1 ) + J , ( 2 )/T|( 2 ) 
are obtained from the orthotropic material constants 


(54) 


s = 


Cl 122 + 


VC1111C2222 


-1 


(55) 


T] = 


C1212 (\/Ci 111C2222-C1122) 
C2222 (2C1212 + Ci 122 + \/Ci 1 1 1C2222) 


1 1/2 


fCn 

V C 22 
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C2222 


where Qju is the stiffness matrix that relates the stresses, <7, j, to strains, £*/ 

3 2 

&ij = ^ ^ Ct /fc/Efcf 
*=1/=1 


(56) 


(57) 


(58) 


By substituting Eq.49-50 into Eq.26-30, the simple relationship between the stress intensity 
factors and the energy release rates are obtained: 


011*121+012 022021 - 
L 022021 — O 12 011021 + 


« 2 n] / M? U 

o ?2 J 1 Pi*?; J l 


021 A,] G/ — o 1 2 A. 2 G/ 
021 A-i G// + 012 X 2 G/ 


-// J 


(59) 


where 


8 cosh 2 (tte) (l +4e 2 ) P 2 
sO) — 5 ( 2 ) 


(60) 


, 16 cosh 2 (7te) (l +4 e 2 ) P 3 /JS1<k 

2_ (sW-sW)( p,+p 2 ) ( } 

The square of the stress intensity factors, Kf and Kf { , can be easily solved from the linear 
equation in Eq.59. The signs for the stress intensity factors, Ki and Ku, can be ascertained from 
the crack opening displacements. Once again, for a homogeneous orthotropic continuum, the 
relationship between the stress intensity factors and energy release rates are reduced to 
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Kj = — G/, Kjj = 2ri G// 

5p 


( 62 ) 


3.4 Interfacial crack in an anisotropic bimaterial continuum 

For a general anisotropic bimaterial continuum, the asymptotic stress field along the interface at a 
distance r ahead of a crack tip is given by 


1 


and the crack opening displacements at a distance r behind the crack tip are 


(63) 


5(r) = ^2r/nDY 


(64) 


[(1 +2/e)cosh (xce) 

where T = {cr 2 i, a 22 , a 23 } r , 8 = {5i, $ 2 , 5 3 } r and { K n , K, Km} T . The bimaterial matrices 
D and W can be written in terms of the anisotropic matrices L and S (Eqs. 5 and 5) 


D = L (1)_1 +L (2) ~' 


(65) 


W = s( 1 )l (1)_1 -S( 2) L (2 ) _1 

and the bimaterial function Y defined 6 in can be rewritten as 


( 66 ) 


Y (r* £ ) = v(^ + sin(elnr)V( 2 ^+cos(elnr)V( 3 ^ 


DY 


r* e 


(1 + 2/e)cosh (7te)J 


= + sin (elnr) U ^ + cos(e!nr) U ^ 




in which the matrices 


(67) 

( 68 ) 

(69) 


V W=I+^(D-'W) 2 , VW = -\D- X W, V( 3 ) = (D~'W) 

(JW =D+±WD~ l W, U® = 
t/( 3 ) = 


(l-j-4e 2 ) cosh(yce)p [ 


(l+4e 2 ) cosh(rce)P 


2eW-^WD~ l W 


I 

W+fWD-'W 


(70) 


24 


A procedure similar in the case of an isotropic bimaterial continuum can also be applied in the 
case of an anisotropic bimaterial continuum. For a general anisotropic bimaterial continuum, three 
additional energy release rates are required, 


Gm = ttSo a 23 ( r )^3 (A — r) dr 


2 A 


« e*0) < 

+F 3[y]t 


F 3 ( W {'" ( M 3H - M 3m) + '« («?[}] - *4 m) } 

l/ 21 (4'[L] ““ML]) + ' 22 { U m ~ M 3 2 [Z]) }] 


(71) 


G n- n i = 


T-r [ [C \2 (r) 5 3 (A - r) + a 23 (r) 5i (A - r)]dr 
2A Jo 

= ” 2A K W {** 1 (“3W “ M 3m) + '12 (4pj “ W 3 2 [/]) } 


+F w\ {' 21 ( M 3 H - *4 w) + ' 22 ( M 3W “ u m) } 
+F 3[*1 {' n (“iH" u 3h) +'12 (“iw -*4w) } 

+F m {' 21 (“i w ” 4 m) +' 22 (“![/] “ w) }] 

Gi-m = [a 22 (r) 6 3 (A - r) + a 23 (r) §2 (A - r)] dr 

= [ F 2 [*l {'» Kh ” “Mh) + tn (“3W “ u m ) } 


+F, 


+F, 


Two additional integrals are introduced, 


r %/2 

/ 3 = / cos 2 0(cos0) 2ie dQ i 

Jo 


2 M {^21 ( w iVm] W 3 2) m]) + '22 (“^[Z] *4 w) } 

+ f 3(.i {'ii (*4h -*4h) +' 12 (4m -4 2 !]) } 

3 |j) {' 21 (“2 w - 4 h) +'22 (4w - 4 2 !]) }] 


7t 17 2 3 . 

-7 - - TT& 


100 


10 


(72) 


(73) 


(74) 


(75) 


/•ti/ 2 1 /jr IQ \ 

*-j£ cos2e(sine) rfe “3SRii)(4-To i£ ) < 76) 

The stress intensity factors can be related to the energy release rates through the matrix A; 


3 3 




n= 1 m— 1 


(77) 
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where 


= 

^2 = 


1 

2kA 

1 

2tcA 
1 


ft /?*- 1 

f A M — t 

J sin(elnr)y dr 


= — [cos (e In A) Im /4 + sin (e In A) Re/ 4 ] 


= - [cos (£ In A) Re/4 — sin (e In A) Im/4] 

7t 

54 = 

= — [cos(elnA)Im/3 -f sin (e In A) Re/3] 

K 

= 2maIo sin ( elnr ) sin ( eln ( A “ r )) \j 

= — [— cos (2eln A) Re/i + sin (2eln A) Im/j + Re/2] 
2 71 

1 /a — r 

^6 = 2kAJo cos ( e ln r ) sin ( e ln C A — Y — ~ 

1 

[cos (2eln A) Im/i + sin (2elnA) Re/j +I111/2] 


2n 


= 2SA / A C ° S (e ln (A - ? ) 

= — [cos (e In A) Re /3 — sin (e ln A) Im/ 3 ] 

K 

1 I A — T 

^8 = 2 ^aJ 0 sin (elnr) cos (eln(A — r)) y - dr 


2n 


[cos (2eln A) Im/i + sin (2elnA) Re/i — Im/ 2 ] 


* - u 


cos (e ln r) cos (e ln (A 


~r))^r 


(78) 


= — [cos (2e ln A) Re/i — sin (2e ln A) Im/i + Re/ 2 ] 

2n 

The matrix B which relates the stress intensity factors k with the energy release rates g using 
the matrix A. 

Bk = g (79) 


where 
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k = {Kj, Kj h Kj Ih K„K nh KjKf/i, K,K n } (80) 

g = {G/, Gui Gui , Gii-iih Gi-ui, Gj-u] ( 8 1 ) 



0(2222 

0(1221 

0(3232 

0(2112 

0(1111 

0(3113 

0(2332 

0(1331 

0(3333 

0(2132 + 0(2312 

0(1131 +0(1311 

0(3133+0(3313 

0(2232 + 0(2322 

0(1231 +0(1321 

0(3233 + 0(3323 

0(2122 + 0(2212 

0(1121 +0(1211 

0(3123 + 0(3213 

0(1223 +0(3221 

0(2223 + 0(3222 

0(1113+0(3111 

0(21 13 + 0(31 12 

0(1333 +0(3331 

0(2333 + 0(3332 


OCi 133 + 0(3131 +OCi3i3 + 0(3311 0(2133 + 0(3132 + 0(2313 +0(3312 

0(1233 +0(3231 +0(1323 +0(3321 0(2132 + 0(2312 + 0(2132+0(2312 

0(1123 + 0(3121 +0(1213 + 0(3211 0(2123 + 0(3122 + 0(2213 + 0(3212 

0(1222 + 0(2221 
0(1112 + 0(2111 
0(1332 + 0(2331 

0(1132 + 0(2131 +0(1312+0(2311 
0(1232 + 0(2231 + 0(J322 + 0(2321 
0(1122 + 0(2121 +0(1212 + 0(2211 


(82) 


(83) 


The square of the stress intensity factors, Kj , K] } and Kj n , can be easily solved from the linear 
equation in Eq.79. The signs for the stress intensity factors, K { , Ku and Km, can be ascertained 
from the crack opening displacements. 


3.5 Numerical examples 

Symmetric laminates under tensile loading 

Consider a symmetric laminate under uniaxial tension as shown in Fig. 8 . The edge delamination 
crack front is assumed to be parallel to the *3 axis and growing in the *2 direction along the 
interface between two dissimilar laminae. Therefore, only the cross section along the specimen 
width is modeled using the quasi three-dimensional finite element method. For a generalized plane 
strain condition, the displacement of the “Quasi 3-D” deformation has the form of 
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x 2 



£33 X 3 


Figure 8: Schematic of edge delamination cracks in a symmetric laminate under uniaxial strain 



Figure 9: The mesh for [± 45/0/90]* and [0/±45/90]* laminate with an edge delamination crack 
between the 0 and 90 degree ply 


mi = U\(x\,x 2 ) 

«2 = U 2 (x\ ,x 2 ) (84) 

m 3 = U 2 (x\,x 2 ) +633x3 

Due to the symmetry conditions of the laminate and the uniaxial load, only a quarter of the 

cross section needs to be modeled (with the mesh shown in Fig. 9). 

Two different laminates (made of T300-5208 carbon epoxy) are considered in the present study; 
[0/± 45/90] 5 and [±45/0/90] 5. The onset of delamination crack propagation in both laminates 
can be attributed to the mismatch in the lateral contraction (due to the Poisson’s effect) between the 
90 plies and plies of different angles. The computed stress intensity factors with the ply thickness 
as the characteristic length are tabulated in Tab. 5. This result using the virtual crack closure 
method compares favorably with the result by Chow et al. (1995) using the mutual integral method 
with the particular solution. The difference between these two solutions is less than 2%. 
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Table 5: Stress Intensity Factors for a delamination crack in a laminate under uniaxial tension 


Laminate 

Crack 

between 

Virtual Crack Closure 
Integral Method 

£33 £33 Kr £33^/// 

Kjc Krrc Kmc 

Interactive Integral with 
Particular Solution 

£33 Ei £33 ju tjiim 

Kic Kuc— Kmc 

[±45/0/90] s 
[0/ ±45/90] s 

0/90 

-45/90 

145.9 29.0 -1.4 
80.4 64.9 -6.9 

145.2 29.5 -1.4 
79.5 65.6 -7.0 


Characteristic length = t p i y 

K/c = 0.88 MPa y/m , K IIC = 2.9 MPa y/m , and K mc = 2.3 MPa^/m. 
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